First Passage Time Densities in Resonate-and-Fire Models 
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Motivated by the dynamics of resonant neurons we discuss the properties of the first passage 
time (FPT) densities for nonmarkovian differentiable random processes. We start from an exact 
expression for the FPT density in terms of an infinite series of integrals over joint densities of level 
crossings, and consider different approximations based on truncation or on approximate summation 
of this series. Thus, the first few terms of the series give good approximations for the FPT density 
on short times. For rapidly decaying correlations the decoupling approximations perform well in the 
whole time domain. 

As an example we consider resonate-and-fire neurons representing stochastic underdamped or 
moderately damped harmonic oscillators driven by white Gaussian or by Ornstein-Uhlenbeck noise. 
We show, that approximations reproduce all qualitatively different structures of the FPT densities: 
from monomodal to multimodal densities with decaying peaks. The approximations work for the 
systems of whatever dimension and are especially effective for the processes with narrow spectral 
density, exactly when markovian approximations fail. 

PACS numbers: 05.40.-a, 02.50.Ey, 87.17.Nn 



I. INTRODUCTION 

The first passage time (FPT) is the time T when a 
stochastic process x(t) leaves an a priory prescribed do- 
main A of its state space for the first time, assumed that 
x{t) has been started at t = from a given initial value 
within A. This concept was originally introduced by E. 
Schrodinger when discussing behavior of Brownian par- 
ticles in external fields Q]. A large variety of problems 
ranging from noise in vacuum tubes, chemical reactions 
and nucleation to stochastic resonance Qjbehavior of 
neurons Q , and risk management in finance [f| can be re- 
duced to FPT problems. In majority of applications the 
attractor of the system's dynamics lies inside A. The 
escape process is characterized by the noise-induced flux 
through the absorbing boundary of A, i.e. by the prob- 
ability density T(T) of the first passage time. 

Approaches to find F(T) are typically based either on 
the Fokker-Planck equation with an absorbing boundary 
or on the renewal theory 0. Despite the long his- 
tory, explicit expressions for the FPT density are known 
only for a few cases. These include overdamped particles 
under the influence of white noise in the force free case, 
under time-independent constant forces and linear forces 
[J IE E3 j as wen as the case °f a constant force un- 
der colored noise [TH ]. Reasonable approximations exist 
for a few nonlinear forces |l2|, [l3|. The FPT densities 
of stationary markovian processes have a very habitual 
form: T(T) goes through a single maximum and then 
it decays either exponentially or as a power law. Many 
neuronal systems do demonstrate such kind of behavior. 
It was shown already in that the interspike interval 
(ISI) histograms obtained experimentally from output of 
some neurons can be reproduced by FPT densities of one 
dimensional diffusion process. 

This kind of description is suitable for overdamped sys- 
tems, where the relaxation time to the attractor t re i is 



much smaller than the typical first passage time. At first 
the local quasiequilibrium is established in the system. 
The escape occurs then from this equilibrium state and 
follows with a constant rate k inversely proportional to 
the mean FPT. This situation is closely related to the 
Kramer's problem considering the quasistationary flux 
over transparent boundary in the low noise limit. The 
problem is independent of the detailed initial state and 
of the time the trajectory has spent inside A. At times 
T exceeding t re i the FPT probability density decays ex- 
ponentially: T(T) ~ exp(— kT). Well known examples 
are chemical systems and nucleation processes, where the 
rates determine the mean velocity of chemical reactions 
or of forming overcritical nuclei [2( • Another example are 
the leaky integrate-and-fire and similar neuronal models, 
where after the reset the corresponding trajectories ap- 
proach quickly the stable rest state [l4[ . 

However, if the time scale separation between the re- 
laxation and escape does not hold, the escape can occur 
before the establishment of the quasiequilibrium and the 
rates are time dependent. The first passage time depends 
sensitively on the initial conditions and the FPT densi- 
ties have a complex shape different from an exponential 
decay. In particular this is the case on short time scales 
T < t re l, which attracted growing interest because of 
recent experiments studying chemical reactions on time 
scales down to femtoseconds. The flux over the bound- 
ary before the establishment of the quasiequilibrium was 
found to grow in a stepwise manner, for an underdamped 
potential system staying initially at the bottom of the 
well [l5|. 

Our work is mainly motivated by dynamics of resonant 
neurons 0, 0, 0] . The voltage variable of such a neu- 
ron exhibits damped subthreshold oscillations around the 
attractive rest state. The characteristic relaxation time 
to the rest state is large compared to the mean ISI. The 
escape of the voltage over the excitation threshold is the 
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beginning of a new spike. After spiking the voltage vari- 
able is reset to a fixed value far from the rest state and 
then it can reach the threshold prior relaxation to the 
rest. The interspike interval densities obtained from the 
output of resonant neurons show a sequence of decaying 
peaks separated by intervals whose length is of the order 
of the period of subthreshold oscillations. In contrast to 
the Kramer's rate theory we stress again the nonstation- 
ary character of this problem due to the reset to sharp 
initial conditions. 

The multimodal ISI probability densities can be re- 
produced in models with different mechanisms of sub- 
threshold resonance: in Hodgkin-Huxley model [l^, in 
excitable FitzHugh-N agum o model with the stable fix 
point bei ng a focus [2(J,[21j or in the region of a canard bi- 
furcation |'2'2tl2.'jl| . It was also shown that two-component 
approximations might well mimic the spiking activity of 
the stochastic Hodgkin-Huxley system [2J, |25( . All these 
models have in common that the multimodal FPT den- 
sity is obtained for stochastic dynamical systems, which 
have at least two dynamical variables, exhibit weakly, 
moderately damped or self-amplifying oscillations and af- 
ter a spike reset to initial values which are not a fixed 
point. In the noise- free situation these systems were de- 
noted resonate-and-fire neurons j2(J . 

In the present work we aim to model excitable be- 
havior with damped subthreshold oscillations. First we 
present the general exact expression for the FPT den- 
sity for stochastic processes with differentiable trajecto- 
ries 0, I27], l28l I23. It results in an infinite series of 
integrals over joint densities of multiple level crossings. 
The later sum stands for a sequential summing of trajec- 
tories excluding all except the ones yielding the first pas- 
sage. Furthermore, we discuss approximations for FPT 
densities, which are based either on truncation of this 
series, or on its approximate summation based on decou- 
pling. We prove the quality of different approximations 
by explicit calculations for an underdamped harmonic 
oscillator driven by white or colored Gaussian noise, rep- 
resenting the stochastic resonate-and-fire neurons. 



II. EXACT EXPRESSION FOR THE FIRST 
PASSAGE TIME DENSITY 

Consider a single random variable x(t), whose t- 
dependence is assumed to be differentiable. The first 
passage problem for x(t) to a boundary Xb is a special 
case of a level crossing problem. 

The general theory of level crossings by a random pro- 
cess was originally put forward by S.O. Rice (2l| . He 
derived an expression for the probability density of re- 
currence of a stationary random process to a given level 
in form of the so called Wiener-Rice series |28[ . The exact 
expression for the first passage time probability density 
Eq.Q is analogous to the Wiener- Rice series and was 
discussed in [2!j, where the main result, our Eq.0, was 
proved. We proceed by giving a much more elementary 



derivation of Eq. Q which serves as the main instrument 
in our further investigations. 

Let us first discuss the probability ni(xb, t\xo, vo)dt 
that a continuous differentiable process x(t) crosses the 
level Xb in a time interval between t and t + dt with pos- 
itive velocity v{t) — x{t) > under initial conditions 
x(Q) — xq,x(0) — vq. Generally the whole set of vari- 
ables resulting from the markovian embedding of x(t) 
should be given at t = 0. For simplicity we consider 
two dimensional dynamics, generalization for the higher 
dimensional systems is obvious. Crossing the level with 
positive velocity will be referred to as an upcrossing in 
what follows. 

If x(t) crosses the barrier within time interval (t, t + dt) 
with velocity v > 0, then the value of coordinate 
at time t should lie in interval Xb — vdt < x(t) < 
Xb- The probability that x(t) is in this interval equals 

J P(x,v,t\xQ,VQ,Q)dx — \v\P(xb, v, t\xQ, t>o, 0)dt. 

Xb—vdt 

Now, the velocity value at the instant of crossing is pos- 
itive but otherwise arbitrary. Thus we obtain the prob- 
ability density of an upcrossing by integration over all 
positive v: 

00 

m(x b ,t\xQ,v ,<$) = J vP(x b ,v,t\x o ,vo,0)dv. (1) 


Eq. can be simply generalized to give the expression 
for the joint probability density of multiple upcrossings. 
The probability n p (xb, t p ; . . . ; x b ,ti\xo, vq, 0)dt p . . . dti, 
that the process x(t) crosses the level Xb in each of p 
time intervals {t\,ti + dti), . . . , (t p , t p + dt p ) is given by: 

00 00 

n p (xb,t p ;...;x b ,ti\x ,v ,0) = J dv p . . . J dv x 



v p . . . vxP(x b , v p , t p ; . . .;a;&,i!i,ti|xo,t>0)0). 

In what follows we omit Xb and initial conditions in 
expressions for the joint densities of upcrossings. The 
transition probability densities are connected with joint 
probability densities according to Bayes' theorem: 

P(x p ,v p ,t p ; . . . ;xb,Vi,ti\xo,v , 0) = 
P2 P +2(x p ,v p ,t p ; . . .;xb,vi,t.xo,vo,0) (3) 
P 2 (x ,v ,0) 

Our aim now is to calculate the first passage time prob- 
ability density J-{T), that is the fraction of all trajectories 
starting from the initial point xo with initial velocity i>o 
which perform the upcrossing of the barrier at time T 
and this upcrossing is the first one. All such trajectories 
are accounted for in probability density «i(T) (see the 
first row in Fig. However, n\(T) also accounts for 
trajectories for which the upcrossing at time T was not 
the first one, i.e. which had another upcrossing at some 
earlier time t\ < T (row 2 in Fig. Such trajectories 
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FIG. 1: Counting crossings. The JVth row corresponds to trajectories with exactly TV upcrossings. The pth column corresponds 
to the pth term of the sum Eq.@. The number on their intersection gives, how many times a trajectory with exactly N 
upcrossings is accounted for in the pth term. Sum of all numbers in every row is exactly zero. 



should not contribute to T(T), therefore we should sub- 
stract them from ni(T). Taking into account that t\ can 
be arbitrary between and T, we get: 



ni(T) - / n 2 (T,t 1 )dt 1 



(4) 



This excludes all trajectories which cross Xb exactly twice 
until T . However, Eq.Q does not fully solve the problem 
since the trajectories crossing Xb three times, i.e at time 
T and at two earlier moments ti < T,i — 1,2 (row 3 
in Fig. nj, are not accounted for correctly. Each such 
trajectory is counted once in n\(T). The second term 
in Eq.(@J accounts for the pairs of upcrossings at T and 
at some U < T . Each trajectory with two additional 
upcrossings at i, < T, i = 1, 2 is therefore subtracted 
twice in J Q n 2 (T,ti)dti. Such trajectories should not 
contribute to T{T): in Eq. (QJ we subtracted too much 
and have to add the amount of trajectories with three 
upcrossings at times < U < T, i = 1, 2 and T again: 



T T 



ni(T)- j n 2 {T,t 1 )dt 1 - 
o 



n 3 (T,t 2 ,ti)dtidt 2 . (5) 



o o 



The factor 1/2! in the last term accounts for the 
number of permutations of variables ti. Generally, 
if a trajectory crosses the level at time T and at 
N earlier times ti < T,i = 1,...,N, then in 
■h J ■ ■ ■ J n p +i(T, t p , . . . , t\)dti . . . dt p it is accounted for 



exactly C P N times (C V N stands for the number of com- 



ix' 



0. 



binations). Note, that J2 {-l) p C p N = (1-1 

Thus in the alternating sum of the kind Eqs.(QJ, JSJ) con- 
taining N + 1 terms, all trajectories crossing Xb at time 
T and having i = 1, 2, . . . , N additional upcrossings are 
excluded, however the ones with the larger number of up- 
crossings are not accounted for correctly. Extending the 
sum to infinity we exclude all superfluous trajectories, 
and only trajectories, for which the upcrossing at time 
T was the first one, remain. Thus the expression for the 
first passage time probability density reads: 



HT) = E 



p=0 



(-1)* 



n p +i(T, tp, . . . , ti)dt p . . . dt\ 



Explicitly expressing the joint densities of upcrossings 
using Eqs.©,©, we get 



T(T) 



- OO 

i V 



(-1)* 



P 2 (xo,v ,0)^ pi 

P-u o 



dtp . . . dti 



OO OO 



(7) 



dvdvi . . . dvpWi . . . v p 



P2 P +4(xb,v,T; x b , v p ,tp] . . . ; x b , vi,ti;x , v Q , 0). 

Eq.© connects F{T), i.e. the solution of the FPT 
problem with absorbing boundary at Xb, with all joint 
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densities of upcrossings for the unbounded process. To 
obtain n p (t p , . . . ,t±) we consider trajectories, which are 
not absorbed at Xb, but can return after an uncrossing 
and then cross Xb again and again. The right combination 
of all these densities of multiple level crossings results in 
the probability density for the first upcrossing. 

An alternative to direct summation of infinite series 
Eq.© is based on an analogue to a cumulant expan- 
sion. The times, when the random process x(t) performs 
upcrossings of Xh form a point process, or a system of 
random points |7], lid ] . The functions 



ni(ii 



n 2 (t2,h 



n 3 (t 3 ,t 2 ,ti 



(8) 



are the distribution functions of the point process. Since 
x(t) has finite velocity, the interval between two upcross- 
ings can not be arbitrary small and so n p (t p , ... ,ti) van- 
ishes if two of its arguments coincide. Such random point 
process is called a system of nonapproaching points |10| . 
In this context J~{T) is interpreted as the waiting-time 
density of the point process. The last is the probability 
density for the time T when the first event occur. 

The system of random points is completely character- 
ized by its cumulant functions 



gi (*i), 52(^2,^1), 53(^3,^2,^1), 



(9) 



Choose an arbitrary natural number r, and then fix r ar- 
bitrary numbers z\, . . . , z r and r positive times t\, . . . , t r . 
The cumulant functions arc then defined by the relation 



1 -\- ^ ^ - ^ ^ n p (t ai . . . , t^Za . . . z u 



-1 a,...,u)=l 



exp 



/ - 1 ' \ 

I / J 1 / t 9p (tat , • • • , ) Zct ■ ■ ■ %w I 
\p=l P' a,...,uj=l / 



(10) 



We obtain the explicit relations between cumulant and 
distribution functions of the point process if we differ- 
entiate both sides of Ea. (|lU|l over all Z{ and then set 
Zi = 0, (i = l,...,r), i.e. if we apply the operator 
d r /(dzi . . . dz r )\ z _ _ z _ Q . Doing so sequentially for 
r = 1, 2, 3, . . . we get 



9i(ti) = niih), 

g 2 (t 2 ,h) = n 2 {t 2 .t 1 ) - ni(ti)ni(f 2 ), 
9z{h, t 2l ti) = n 3 (i 3 ,t 2 ,ti) - 3{ni(ti)n 2 (t 3 ,t 2 )} s 
+ 2n 1 (t 1 )n 1 (t 2 )n 1 (t 3 ), 



(11) 



Here {. . .} s denotes the operation of symmetrization of 
the expression in the brackets with respect to all permu- 
tations of its arguments. The coefficients in these forms 
are the same as in relations between the moments and 
the cumulants of a random variable. 

The relation Ea. docs not change its form if we 
choose different times t\, . . . ,t r , different values zi, . . . ,z r 
or change the number r. Thus extending r to infinity, 



allowing t to take all possible values between and T, 
and choosing Z\ =...= z r = —1 we get from Eq. l|lU|) 



P =i 



(-if 



T T 



n p (t p , . . . , ti)dt p . . . dti 







(12) 



exp J2 



T T 



g p {t p , . . . , ti)dt p . . . dt\ 







It is easy to verify, that the derivative (—d/dT) of the 
expression on the left hand side of Eq. i|12fl is exactly the 
expression on the right hand side of Eq. © . Thus differ- 
entiating the right hand side of Eq. I|12() over T we obtain 
the expression for the waiting-time density T(T) through 
the cumulant functions of the point process: 



T(T) = S'{T)e- s{T) 



(13) 



with 



5 ( T ) = -E^ ! j ■■■j'M',- 
P— 1 n n 



, t\)dt p . . . dt\. 



(14) 

Eas.(P|l- H12|) are general expressions which hold for sys- 
tems of random points, defined by distribution functions 
Eq.© of any kind. So are also Eqs.® and (fT5]b. ((Tlj) . 
which give the waiting-time density for arbitrary point 
process. In particular, for the random points being the 
times, when a differcntiable random process crosses the 
level Xb, the distribution functions n p (t p , . . . ,t%) are given 
by the joint densities of upcrossings Eq. © . Then Eqs. © 
and i|13|) , (|14fl together with Eq. @ express the first pas- 
sage time density for this differentiable random process. 
The function S"(T) can be interpreted as a the time- 
dependent escape rate. 

These are the exact results for the FPT probability 
density of any continuous differentiable random process. 
Though these results were employed for mathematical 
proofs, up to our knowledge these infinite series of mul- 
tiple integrals was never used for explicit calculations. 
We proceed to show, that Eqs.© and ljT3|) . (|Tlfl can be 
a starting point for several approximations. As often in 
the case of infinite series, the useful approximations can 
be based either on the truncation of the series after sev- 
eral first terms calculated exactly, or by approximation 
of the higher order terms through the lower order ones 
what might lead to a closed analytical form. Truncation 
approximations for Eq.JfjJ are not normalized, hold only 
on short time scales, and diverge at longer times (due 
to the miscount of trajectories with several upcrossings). 
The approximations of the second type are based on a 
subsummation in Ea. H14|) for S(T). They are normalized 
and can be used in the whole time domain. Note, only 
approximations guaranteeing positive rates S'(T) are rea- 
sonable. Thus the set of possible approximations for the 
series Ea. (|14|l is rather restricted. 
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III. NOISY DRIVEN HARMONIC 
OSCILLATOR: RESONATE-AND-FIRE 



The model we have in mind is the resonate-and-fire 
model of a neuron [2(|. This is the least complicated 
model accounting for the resonant properties of neurons 
in terms of an equivalent RLC circuit [l7j . In this 
way it is directly related to the leaky integrate-and-fire 
model also based on the electrical analogy. Alternatively 
the model can be interpreted as a systematic and lin- 
earized reduction of Hodgkin-Huxley type dynamics |3f)j | . 
For the sake of simplicity we neglect the absolute refrac- 
tory time. Under this assumption the model is equiv- 
alent to the underdamped harmonic oscillator with the 
threshold and reset, what makes the results applicable in 
many other domains of science. We change by time scale 
and variable transformations to dimensionless parame- 
ters and variables. The dynamics of the voltage variable 
x{t) is given by: 



v = —71; — u) x + rj(t). 



(15) 



We fix the frequency u>q = 1, choose initial conditions 
for x and its velocity t) to be io — — 1j w o — and 
set the threshold at Xb = 1. In the present paper we 
consider two types of noisy drive: (i) the white noise 
■q{t) = \/2Dt;(t), and (ii) the Ornstein-Uhlenbeck noise 
fj = -T- 1 !] + VzDt^ 1 ^), with £(t) being the white 
Gaussian noise of intensity 1. First we concentrate on 
the case (i) of white noise driving. 

Because of the linearity of the system Ea. (|15fl all joint 
probability densities are Gaussian and have the form : 



Pn(Q) 



(2tt)«/ Vdet C n 



exp 



QC~ X Q 



(16) 



Here Q = (qi(ti), . . . , q n (t n )) is an n-dimensional vector, 
whose zth component is the value of coordinate x(ti) or 
of velocity v(ti) at the moment U. C n is symmetric nxn 
correlation matrix. Its elements are correlation functions 
between corresponding components of vector Q : Cij — 
cji = (qi{ti)qj(tj)). 

Correlation functions for the system Eq. 1|15|) are easily 
obtained using Fourier transform and Wiener-Khinchin 
theorem |6j. For the case of white noise driving 
and in an underdamped regime (7 < 2ujq) r xx (t) = 
(x{t')x(t' + t)) = -jpre-i* (^sin(fit) + cos(fii)) 5 with 



Q 



21 

4 



In overdamped case (7 > 2uo) the 
expression for r xx (t) is the same, except the trigonomet- 
ric functions are replaced with hyperbolic ones. Further, 
r xv{t) = r' xx (t) and r vv (t) = ~r^ x (t). 

Then n\(T) is obtained from Eq.QJ in closed analytical 



form: 

m(T) : 
exp 

Here 



27r/i22 V det C4 
- 9 £ /'<.•'/<'// 



exp 



r cr 2 x vl +a 2 v xl 



2a x (T v 
1 — y/nae" erfc(a} 



(17) 



(J>ji 



are elements of the inverse correlation 



matrix (C4) , qi are components of the vector Q = 
(x(T), v(T), £o(0), ^o(O))- The dispersions of x and v are 
o 2 x = r xx (0) = D/ 7 wg, a 2 v = r vv (0) = D/j. We have 



introduced a= ^] \inq% /V^22- Finally, erfc(a;) is a 

W 2 / 

complementary error function. 

For the joint densities of multiple upcrossings 
n p (t p , . . . , ti) no closed expressions can be obtained. We 
evaluate the integral over v\ in Eq.0 analytically and 
then perform numerical integration of the resulting ex- 
pression over «2, . . . , v p to obtain n p (t p , . . . , ti). The in- 
tegrals over time in the expressions for T(T) are also 
evaluated numerically. 



IV. TRUNCATION APPROXIMATIONS 

The first passage time density T(T) for the harmonic 
oscillator driven by white Gaussian noise obtained from 
simulations is depicted with black crosses in FigEJ Pa- 
rameters are chosen to be 7 = 0.01, D = 0.02. In this 
case of very small friction the correlation fu nctions of 
the process oscillate with period T p — 2n/ y/u;^ — 7 2 /4 
and decay slowly within the relaxation time t re i = 2/ 7 . 
A typical trajectory is smooth and shows almost regular 
oscillations with fluctuating phase and amplitude. The 
probability to reach Xb is higher in the maxima of the sub- 
threshold oscillations. The initial phase of these oscilla- 
tions is fixed by sharp initial conditions. Thus on shorter 
time scales T(T) shows the multiple peaks following with 
the frequency of damped oscillations = \/uJq — 7 2 /4. 
On long times T ^> t re i the quasiequilibrium establishes 
and FPT density decays exponentially. The number of 
visible peaks depends on the relation between t re i and the 
period of oscillations T p and is given by the number of 
periods elapsing before the quasiequilibrium is achieved. 
For parameter values as in Fig|2trei = 200, what corre- 
sponds to about 30 periods T p = 6.28. 

Let us consider truncation approximations for the se- 
ries Eq.©. The first approximation is given by the 
first term n\{T), the second approximation by 2 terms, 
Eq.(gJ, and the third by 3 terms, Eq.©. The higher 
order approximations entail the numerical estimation of 
highdimensional integrals, which at some stage leads to a 
computational effort larger than the one necessary for a 
direct simulation. Therefore we restrict ourselves to the 
1-, 2- and 3- terms approximations. 
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cess, also known as the Rice frequency [2g ■ The general 
expression for rig reads 



n 



1 

2^ 



,(0) 



,(0) 



1/2 



= -x^/2r :EX (0) 



(18) 



In our case of a harmonic oscillator driven by white 
noise n = (wq/2-k) exp (— jxIujq/2D). In the stationary 
regime the mean interval between two consecutive up- 
crossing Tfj is given by the inverse of the Rice frequency 
Tr = 1/riQ. For the chosen parameter values Tr = 8.06. 

The second approximation (gray line in Fig^b)) re- 
produces almost exactly the first two peaks of FPT den- 
sity. Then it becomes negative, because in Eq.|@} tra- 
jectories performing 2 and more superfluous upcrossings 
are subtracted too many times. Moreover the second 
approximation tends to minus infinity for T — > oo. The 
third approximation reproduces well the three first peaks 
of F{T), and then diverges tending to plus infinity. 

Note that the mean first passage time obtained nu- 
merically equals 14.6 for these parameter values, and 
the median of the distribution lies by 3.2. Thus the 
first three approximations reproduce the most part of 
the FPT probability density. 
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FIG. 2: FPT density for harmonic oscillator driven by white 
Gaussian noise, u>o = 1, 7 = 0.01, D — 0.02, xo = — 1, Vo = 
0,xt — 1. Simulation results shown with black crosses and 
truncation approximations with gray line: (a) 1 term Eq.Q, 
(b) 2 terms Eq.lgJ and (c) 3 terms Eq.©. The mean FPT 
obtained from simulations equals 14.6, the median of distri- 
bution lies by 3.2. Note the logarithmic scale. 



The result of the 1 term approximation is shown in 
Figl^a) with a gray line. The first peak of the FPT den- 
sity is reproduced almost exactly. All further peaks are 
overestimated, because all trajectories performing mul- 
tiple upcrossings of Xf, are included. On long times the 
process becomes stationary and the first approximation 
tends to a constant value limr^ oc ni(T) = uq. This is 
the mean frequency of upcrossings for a stationary pro- 



FIG. 3: FPT density for harmonic oscillator driven by white 
Gaussian noise, ljq = 1, 7 = 0.8, D = 0.44, xo = — 1, t>o = 
0,xt = 1. Simulation results shown with gray crosses; 1 term 
approximation Eq.Q with a dot line, 2 terms approximation 
Eq.JlJ with a dashed line, and 3 terms approximation Eq.JSJ 
with a solid line. The mean FPT obtained from simulations 
equals 13.1, the median of distribution lies by 9.1. Note the 
logarithmic scale. 

The behavior of F{T) for the harmonic oscillator with 
higher damping 7 = 0.8, stronger noise intensity D — 
0.44, and other parameters as in FigEl is presented in 
FigED For these parameter values the relaxation time 
t re i = 2.5 is less than the period T p = 6.86. Therefore 
the FPT density is practically monomodal with a sin- 
gle maximum and a small shoulder separating it from 
the exponential tail. The numerically obtained F(T) is 
shown with gray crosses, the 1 term truncation with a 
dot line, the 2 terms truncation with a dashed line, and 
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the 3 terms truncation with a solid line. The truncation 
approximations reproduce again the most part of the dis- 
tribution: the mean FPT equals 13.1 and the median lies 
by 9.1. 
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FIG. 4: The same as in Fig|5J however the truncations are cor- 
rectly normalized by Eg. 1191 . Simulation results shown with 
gray crosses; the truncations with solid lines. The curves for 
normalized 1,2, and 3 terms truncations are vertically shifted 
by 0.6 for the sake of clarity. 



Thus the truncation approximations reproduce the 
FPT density on short times but are not normalized and 
diverge on large times. However, one can force the nor- 
malization in truncations as follows 



\Ftrunc{T)\ 6 ( 1 - \Tu 



(19) 



Here TtrunciT) is a truncated series for the FPT density, 
it is given by Eqs.Ql,Q), and (JSJ for 1,2, and 3 terms 
truncations respectively. O(i) is the Heaviside step func- 
tion: the expression Eq. (|19J) turns to zero as soon as the 
integral over the absolute value of ^ r t rimc (T) exceeds 1. 

In FigQ] we show the FPT density for the same pa- 
rameter values as in Fig[3 Simulation results are shown 
with crosses, and approximations obtained from Ea. H19(l 
with 1,2, and 3 terms by solid lines. The probability dis- 
tributed over the long exponential tail in the real FPT 
density, is concentrated in the positive artefact posed on 
intermediate times in the normalized truncations. Hence 
the mean FPT computed from such approximations is 
always strongly underestimated. 

The more terms are included, the more precise trun- 
cations become. However one has to confine oneself to 
a few terms, since the calculation of higher order terms 
implies the computation of multiple integrals and is not 
more effective than simulations. Therefore the trunca- 
tion approximations are good, when the most part of the 
FPT probability is concentrated in the first few peaks, 
i.e. when the barrier value is low or the noise is strong. 



V. DECOUPLING APPROXIMATIONS 

Decoupling approximations for Eq.® or Ea. (|14|) are 
based on approximate expressions of the higher order 
terms through the lower order ones, what may lead to 
a closed analytical form. Thereby infinitely many ap- 
proximate terms are included. 

The simplest way to obtain such an approximation is to 
neglect all correlations between upcrossings. This means 
to neglect all terms in Ea. (|14|l except for the first one, 
and leads to 



S{T) = J m{t)dt. 



(20) 



Equivalently, neglecting all correlations corresponds 
to the factorization of n p +i(T, t p , . . . , ti) into a prod- 
uct of one-point densities ni(T)ni(t p ) . . . ni(ti) in 
Eq.©. Then the series, Eq.© sums up into 

T{T) w n\{T) exp J Q T m(t)dt \ , which is equivalent to 

Eas. l|13|l . (|20|) . This approximation will be refereed to as 
the Hertz approximation since the form of J^(T) resem- 
bles the Hertz distribution |3l| . It is an approximation 
of first order, since it takes the first term of the series 
exactly, and all other terms are approximated through 
this first one. 

The second order approximation should therefore ac- 
count the first and the second terms exactly and approx- 
imate all higher terms through these two. The general 
form of F(T) in terms of the cumulant functions Eq. Ijl3(l 
ensures the right normalization, irrespective of the way 
how S(T) is approximated. However the simple trunca- 
tion of the series Eq. Ijl4|l after the second term does not 
ensure the positive escape rate S'(T). 

The second order approximation guaranteeing S'(T) > 
was proposed by Stratonovich in the context of peak du- 
ration 10]. The first and the second cumulant functions 
are taken exactly, and the higher ones are approximated 
by the combinations of these two: 



g p {t p , w(-l)f- 1 (p - l)!ni(tp) . . . m(ti) 

{R(t 1 ,t 2 )R(t 1 ,t 3 )...R(t 1 ,t p )} s 



(21) 



Here is again the operation of symmctrization. 

R(ti,tj) is the correlation coefficient of upcrossings 



R\ti , tj 



n 2 (tj, tj) 
ni(ti)ni(tj) ' 



(22) 



Note, that R(ti,ti) = 1 and R(ti,tj) — » for large values 
of \U -tj\. 

The approximation of the cumulant functions in form 
Eq. (|21|) can be motivated by the following argument. 
Consider Ea. (|10fl with r = 1. Recall that the joint 
densities of upcrossings vanish for coinciding arguments: 
n p (ti, . . . , ti) = 0. Thus it follows from Eq.^TUJl: 

oo 1 

ln(l + m(fi>i) = ^2 —dpiti, ■ ■ ■ M)z{- 
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FIG. 5: FPT probability density for harmonic oscillator 
driven by Gaussian white noise. Simulation results are shown 
with a gray line, Hertz approximation with a black dashed 
line, and Stratonovich approximation with a black solid line. 
Note the logarithmic scale in T. The insets show the same 
curves on the logarithmic scale in T{T). The parameters are: 
uo = l,xo = -l,«o = 0, Xb = 1, a) 7 = 0.8, Z? = 0.1, t re ; = 
2.5, T fl = 343, b) 7= 0.8, D = 0.44, t rei = 2.5, T fl = 15.6, c) 
7 = 0.08, D = 0.01, t rei = 25, T R = 343 



gives a correction to it, when the arguments differ. 

Substitution of Eq. (J2TJ into Eq. fFQ delivers then the 
Stratonovich approximation for T{T) in form Eq. (|13|l . 
now with the time-dependent escape rate being 



^ ' -J »i - 
" 



In 



1 - f Q R^t^mit^dt' 



/ R(t,t') ni (t')dt' 



dt. (24) 



Let us now discuss the domains of applicability for 
these approximations. The Hertz approximation Ea. l|2U|) 
holds if all correlations decay considerably within the typ- 
ical time interval between upcrossings Tr. The decay of 
correlations is described by the relaxation time t re i =2/7 
of the process. Therefore, the Hertz approximation holds 
for t re i < T R - 

The Stratonovich approximation is applicable when 
the argument of the logarithm in Eq.J2U) is positive, 
1 - fT\n->(t') + n 1 (t)- 1 n 2 (t',t)}dt' > 0. Using the fact 



that n 2 (y, t)/ni(t) tends to ni(i') for \t — t'\ > t rei and 
tends to zero for \t — t'\ — > we get as a rough estimate 
for the validity region of Eq. 124(1 t re i < Tr. 
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The above expression should hold for arbitrary Z\. 
Therefore expanding the logarithm in series, and equat- 
ing the coefficients by the same powers of Z\ on the both 
sides, one obtains the identity 

g p (t 1 ,...,t 1 ) = (-iy- 1 (p-l)\n p 1 (t 1 ). (23) 

Ea. (|23ll is exact for all arguments coinciding. Eq. l21|) 



FIG. 6: Same as in Fig. |K| however for the case of stronger 
friction. The parameters are a) 7 = 3.0, D — 0.5, u>o/"f = 0.33, 
b) 7 = 10.0, D = 5.5,uo/7 = 0.1 and other parameters as in 
Fig. 

Let us now turn to the results for the harmonic oscil- 
lator Eo. 1(15(1 with white noise driving. In Figs 03 and 
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the FPT probability density obtained from simulations 
is depicted with a gray line, the Hertz approximation 
Ed. 1)21) fl with a black dashed line, and the Stratonovich 
approximation Ea. (|24l) with a black solid line. 

In Fig[S{a) the parameters are chosen to be 7 = 
0.8,1? = 0.1, corresponding to moderate friction and 
moderate noise intensity. For given parameter values 
t re i = 2.5 and Tr = 343, so that t re i <C Tr, both 
Hertz and Stratonovich approximations hold and repro- 
duce well the FPT density in the whole time domain. 

In the case of moderate friction and stronger noise 
the upcrossings become more frequent and Tr decreases. 
The FPT changes its form to practically monomodal. An 
example is given in Fig^b) with 7 = 0.8, D = 0.44 
which correspond to t re i = 2.5 and Tr — 15.6. The 
Stratonovich approximation complies very well with sim- 
ulations, while the Hertz approximation fails to repro- 
duce the details of the distribution: It underestimates 
T{T) on short times, and shows slower exponential de- 
cay in the tail than the one observed in simulations (see 
the inset). 

Finally, for small friction and week noise the upcross- 
ings are rare, but the relaxation time is large. The FPT 
probability density exhibits multiple decaying peaks. In 
FigG2 c ) 7 = 0.08, D = 0.01 corresponding to t re i = 
25, Tr — 343. Again, the Stratonovich approximation 
performs well, while the Hertz approximation underesti- 
mates the first peak, overestimates all further peaks and 
decays in the tail faster than the simulated FPT density. 

For the large T, T(T) decays exponentially, T{T) oc 
exp(— kT). The decrement of this decay is obtained 
from long time asymptotic: kT — limj^oo S(T). Thus, 
in the Hertz approximation Ea. (|20|l one gets k,h — 
lim T ^ 00 (l/T) f^m(t)dt = n T/T = n . The behav- 
ior in the Stratonovich approximation Ea. H24|) is deter- 
mined by lim^f/^oo J Q R(t, t')n\{t')dt' « n T cor with r cor 
given by r cor = lim f 00 R(t,t')dt'. Note that r cor is not 

necessary positive because of oscillating correlation coef- 
ficient. Inserting this expression into Ea. (|24|l and ex- 
panding the logarithm up to the second term we get 
K s = "o(l + \riQT cor ) providing the second order cor- 
rection to kh- The value of r cor for the parameter 
set as in FigEfa) is r cor = —2.4, for parameters as in 
Fig'EIb) r cor = 5.09, and for parameters as in Fig. \^c) 
T CO r = —431.99. The long time asymptotic obtained with 
these r cor values reproduce fairly well the decay patterns 
found numerically. 

In the overdamped regime (7 > 2wo) the condition 
t re i < Tr is always fulfilled. Nevertheless the validity 
region of our approximations is limited. With increas- 
ing friction the process x(t) approaches to the markovian 
one (it is markovian in the overdamped limit luq/^/ <ti 1). 
For such processes the pattern of upcrossings is not 
homogeneous, but shows rather well separated clusters 
of upcrossings [Tfl j. Essentially in the markovian limit 
the property, that upcrossings form a system of non- 
approaching random points is violated. The upcross- 



ings within a single cluster are not independent even 
if their mean density no is low, so that the quality of 
approximations decreases. This fact is illustrated in 
FigOD I n the overdamped regime the correlation func- 
tions decay monotonously, and the FPT densities are 
always monomodal. The parameters in FigUJa) are 
7 = 3.0, D = 5.5, so that w / 7 = 0.33. Eqs.CSJ,^ 
continue to give a good approximation for !F(T), while 
the Hertz approximation becomes inaccurate. Further 
increase in friction, for example 7 = 10.0, D = 5.5 as in 
FiglUb) corresponding to wq/7 = 0.1, makes the pro- 
cess to approach the Markovian limit. The Stratonovich 
approximation starts to be inaccurate, and the Hertz ap- 
proximation fails. 



VI. TRUNCATION VERSUS DECOUPLING 
APPROXIMATIONS. 

In the two previous sections we have seen, that trunca- 
tion approximations reproduce the FPT density on short 
time scales. In contrast, the decoupling approximations 
reproduce FPT densities in the whole time domain and 
posses the right normalization. At the first glance, it 
may seem that the decoupling approximations excel the 
direct truncations and should be preferably used in ap- 
plications. However, it depends on the problem one has 
to solve, and sometimes the truncations turn out to be 
useful. 



0.25 




FIG. 7: FPT density for harmonic oscillator driven by white 
Gaussian noise, u)o = 1,7 = 0.08, D = 0.01, xa = — l,«o = 
0, xt — 1. Simulation results are shown with crosses, the 3 
terms approximation with black line, and the Stratonovich 
approximation with gray line. The inset shows the magnifi- 
cation of a part of the curve. 

One such situation was already mentioned. If the noise 
intensity is high or the barrier value is low, then the up- 
crossings occur frequently, and the decoupling approxi- 
mations can not be applied. For example, for parameter 
values as in Fig|3 the relaxation time is t re i = 200, and 
the mean interval between upcrossings Tr = 8.07. Thus 



10 



trei > Tr, and both Hertz and Stratonovich approxima- 
tions fail. Nevertheless, the most part of the FPT density 
is concentrated in the first few peaks in this case, and 
is well reproduced by the truncation approximations, as 
shown in Figs and 0] 

One can also be interested in a very accurate approx- 
imation for the FPT density on short times. The trun- 
cations deliver better results on short times than the de- 
coupling approximations. For example, in Fig[7|the sim- 
ulation results are compared with the 3 terms truncation 
and the Stratonovich approximation for the same param- 
eter values as in FiglSJc). Both approximation reproduce 
very accurate the first peak in the FPT density. How- 
ever, the 3 terms truncation is much more accurate in 
estimation of the second and third peaks (see the inset 
in FigEJ. 

This can be easily understood. The Stratonovich 
approximation takes ni(ii) and na^ai^i) exactly, and 
approximates all higher order densities through these 
two. On times, when the first peak occurs, nzfoiti) 
is neg ligibly small. Then from Eas.lfST j) . ^ we ob- 
tain g p (tp, . . . ,h) w (-l)?" 1 ^ - l)\m(t p ) . . .ni(*i). 
Substitution of this expression into Ea.Q10[l gives then 
n p (t p , . . . , t\) w for p > 1. Thus on these times the 
Stratonovich approximation just coincides with the 1 
term truncation, and so reproduces the first peak very 
accurately. On times, when the second peak in the FPT 
density occurs, ^2(^2,^1) is significantly different from 
zero. Hence all approximated n p (t p , . . . ,t\),p > 2 turn 
out to be nonvanishing as well, while the real values 
for these functions are negligibly small on these times. 
Thus the accuracy of the Stratonovich approximation de- 
creases in the second peak. From analogous reasoning it 
becomes clear, that the Hertz approximation is already 
inaccurate in estimation of the first peak. 



the process tends to the white noise of intensity 2D. 
The correlation functions r xx (t),r xy (t),r yx (t),r xv (t) = 

r xx(t)> r vv{t) — ~~ r xx(t)7 r yv(t) — r yx (t)i r vy{t) — ~ r xy{t) 

can be obtained using Fourier transform as it was done 
for the white noise case in Section ITTT1 

Then n\{T) is again obtained analytically and has the 
form given by Eq. (|17fl . Now = fiji are elements of the 
inverse correlation matrix (C5) , qi are components of 
vector Q = (x(T),v(T),x (0),v (0),T]o(0)). The factor a 
is defined in the same way as it was done in Section IlIII 

For simplicity, we assumed sharp initial conditions for 
the noise variable, i.e. 77 is reset after every spike to 
a fixed value t]q. The alternative assumption, that the 
neuron variables x, v are reset to their initial values once 
x(t) reaches the thereshold without resetting n(t), might 
be more realistic ^ n this case all probability den- 

sities should be averaged with respect to the stationary 
density of noise values upon firing. However, as an exam- 
ple, we confine ourselves to consideration of sharp initial 
conditions for the noise: P(r),t = 0) = S(n — rjo). 
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VII. HARMONIC OSCILLATOR DRIVEN BY 
COLORED NOISE 

Expressions Eq.© and Eas.l|13 |) .l|24 fl can be used to 
obtain the FPT density for a random process x(t), if the 
joint probability densities of x and its velocity v Eq.J2J 
exist. Thus it is necessary that the process x(t) is con- 
tinuous and differentiable at any time, but there are no 
further restrictions on dimension and form of the system. 
The truncation and decoupling approximations deliver 
good results for the FPT density in their validity regions 
independently of the character of the noisy drive. In par- 
ticular the case of correlated input signals (colored noise 
driving) is of importance in neuroscience. For example, 
synaptic filtering of the input spike train may lead to an 
exponentially correlated input signal. 

Therefore we consider as another example a resonatc- 
and-fire neuron Ea. (|15f) driven by the Ornstein- 
Uhlenbeck noise. The correlation time of the process 
is r, the variance D/t, and the correlation function 
(r](t)r](t + t')) = {D/t) exp(-i'/r). In the limit r -> 



FIG. 8: Simulation results and the Hertz approximation for 
the FPT probability density for harmonic oscillator driven by 
Ornstein-Uhlenbeck noise with correlation times r = 0.5 (up- 
per curves) and r = 1.0 (lower curves) and other parameters 
as in Figl^Ja). Note the logarithmic scale in T. The inset 
shows the same curves on the logarithmic scale in J-{T). 

In Fig|S| we show simulated FPT probability density 
and the Hertz approximation for the harmonic oscilla- 
tor driven by the Ornstein-Uhlenbeck noise. We choose 
two different values of the correlation time: r = 0.5 and 
t = 1.0. The reset value for the noise is 770 = 0, and 
other parameters are as in Fig[S{a). The noise intensity 
decreases with increasing r, hence the mean FPT growth. 
Though !F(T) preserves its structure as a whole, for larger 
values of r the height of the main peak decreases and the 
weight of the exponential tail growth. For given values 
of the correlation time T R = 829.2 and T R = 4247.1 re- 
spectively, what significantly exceeds t re i. Thus in both 
cases the Hertz approximation is absolutely sufficient. 
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VIII. SUMMARY 

Motivated by studies of resonant neurons, we consider 
the first passage time problem for systems with sub- 
threshold oscillations and nonnegligible relaxation times 
after a reset. The joint densities of multiple upcrossings 
for such a process x(t) can be obtained in the case of dif- 
ferentiable trajectories. The FPT density for x(t) is ex- 
pressed in terms of an infinite series of multiple integrals 
over all joint densities of upcrossings, or equivalently, in 
terms of the cumulant functions. 

We consider two types of approximations for this in- 
finite series. The truncation approximations include the 
first few terms of the series calculated exactly. They re- 
produce well the FPT density on short and intermidi- 
ate times and can be used, when the most part of the 
FPT probability is concentrated in the first few peaks, 
i.e. when the barrier value is low or the noise is strong. 

The decoupling approximations can be derived for the 
case of weakly correlated upcrossings. The higher order 
cumulant functions are expressed through the lower or- 
der ones, and then infinitely many terms sum up to the 
closed expression for !F(T). The Hertz approximation 
(the one neglecting all correlations between upcrossings) 
is absolutely sufficient for the case of moderate friction 
and moderate noise intensity. The Stratonovich approx- 
imation (approximating the higher order cumulant func- 



tions through the first and the second ones) performs 
even better and does not loose the accuracy for high noise 
intensities or in the slightly overdamped regime. 

We illustrate our results by the noise driven harmonic 
oscillator, with the threshold value at Xb and reset to 
sharp initial conditions, i.e. the resonate-and-fire model 
of a neuron. The validity regions of the approxima- 
tions cover all different types of subthreshold dynamics. 
Thus the approximations reproduce all qualitatively dif- 
ferent structures of the FPT densities: from monomodal 
through bimodal to multimodal densities with several de- 
caying peaks. The approximations hold for systems of 
whatever dimension. We illustrate this by the harmonic 
oscillator driven by the Ornstein-Uhlenbeck noise. 

Though we applied the theory to the harmonic oscilla- 
tor (resonate-and-fire model), the linearity of the system 
is in general not required. The joint distributions of x{t) 
and its velocity should exist, i.e. x(t) should be diffcrcn- 
tiable in time. No further restrictions on the form and 
dimension of the system are implied. 
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